The data is simulated as per the writeup here, except that \(\beta_k\), \(k = 1,2,3\) are fixed.
We use Bayesian multivariate generalized linear models with correlated group-specific terms via Stan to model the binary outcomes
In other words, the fitted model was specified as:
\[y_{ijk} = \beta0_{jk} + \beta1_kx_i + \epsilon_{ik}\]
model <- stan_mvmer(
formula = list(
y1bin ~ wealthindex + (1 | years) + (1 | hhid)
, y2bin ~ wealthindex + (1 | years) + (1 | hhid)
, y3bin ~ wealthindex + (1 | years) + (1 | hhid)
)
, data = sim_dflist[[1]]
, family = list(binomial, binomial, binomial)
)
Plots point estimates for each parameter along with credibility intervals. In this case we are plotting the \(50\%\) uncertainty interval (thick horizontal lines) and the \(90\%\) uncertainty interval (thin horizontal lines). In terms of interpretation, the \(50\%\) uncertainty interval identifies where \(50\%\) of the marginal distribution lies for each parameter.